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Abstract 

A new parabolized Navier-Stokes (PNS) algorithm 
has been developed to efficiently compute magne- 
tohydrodynamic (MHD) flows in the low magnetic 
Reynolds number regime. In this regime, the electri- 
cal conductivity is low and the induced magnetic field 
is negligible compared to the applied magnetic field. 
The MHD effects are modeled by introducing source 
terms into the PNS equation which can then be solved 
in a very efficient manner. To account for upstream 
(elliptic) effects, the flowfields are computed using 
multiple streamwise sweeps with an iterated PNS al- 
gorithm. Turbulence has been included by modify- 
ing the Baldwin-Lomax turbulence model to account 
for MHD effects. The new algorithm has been used 
to compute both laminar and turbulent, supersonic, 
MHD flows over flat plates and supersonic viscous 
flows in a rectangular MHD accelerator. The present 
results are in excellent agreement with previous com- 
plete Navier-Stokes calculations. 

Introduction 

Flowfields involving MHD effects have typically been 
computed [1-10] by solving the complete Navier- 
Stokes (N-S) equations for fluid flow in conjunction 
with Maxwell’s equations of electromagnetodynam- 
ics. When chemistry and turbulence effects are also 
included, the computational effort required to solve 
the resulting coupled system of partial differential 
equations is extremely formidable. One possible rem- 
edy to this problem is to use the parabolized Navier- 
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Stokes (PNS) equations in place of the N-S equa- 
tions. The PNS equations can be used to compute 
three-dimensional, supersonic viscous flowfields in a 
very efficient manner [11]. This efficiency is achieved 
because the equations can be solved using a space- 
marching technique as opposed to the time-marching 
technique that is normally employed for the complete 
N-S equations. 

Recently, the present authors have developed a 
PNS code to solve supersonic MHD flowfields in 
the high magnetic Reynolds number regime [12]. 
This code is based on NASA’s upwind PNS (UPS) 
code which was originally developed by Lawrence et 
al. [13]. The UPS code solves the PNS equations us- 
ing a fully conservative, finite-volume approach in a 
general nonorthogonal coordinate system. The UPS 
code has been extended to permit the computation of 
flowfields with strong upstream influences. In regions 
where strong upstream influences are present, the 
governing equations are solved using multiple sweeps. 
As a result of this approach, a complete flowfield can 
be computed more efficiently (in terms of computer 
time and storage) than with a standard N-S solver 
which marches the entire solution in time. Three it- 
erative PNS algorithms (IPNS, TIPNS, and FBIPNS) 
have been developed. The iterated PNS (IPNS) al- 
gorithm [14] can be applied to flows with moder- 
ate upstream influences and small streamwise sepa- 
rated regions. The time iterated PNS (TIPNS) algo- 
rithm [15] can be used to compute flows with strong 
upstream influences including large streamwise sep- 
arated regions. The forward-backward sweeping it- 
erative PNS (FBIPNS) algorithm [16] was recently 
developed to reduce the number of sweeps required 
for convergence. 

The majority of MHD codes that have been de- 
veloped combine the electromagnetodynamic equa- 
tions with the full Navier-Stokes equations result- 
ing in a complex system of eight scalar equations. 
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= V-(V.?)-V-U + E-J 


These codes can theoretically be used for any mag- 
netic Reynolds number which is defined as Re m = 
cTe/ieKooL where <r e is the electrical conductivity, fi e 
is the magnetic permeability, Voo is the freestream 
velocity, and L is the reference length. However, 
it has been shown that as the magnetic Reynolds 
number is reduced, numerical difficulties are often 
encountered [4]. For many aerospace applications 
the electrical conductivity of the fluid is low and 
hence the magnetic Reynolds number is small. In 
these cases, it makes sense to use the low magnetic 
Reynolds number assumption and reduce the com- 
plexity of the governing equations. In this case, the 
MHD effects can be modeled with the introduction 
of source terms into the fluid flow equations. Severed 
investigators [4,8,17-19] have developed N-S codes 
for the low magnetic Reynolds number regime where 
the induced magnetic field is negligible compared to 
the applied magnetic field. 

In the present study, a new PNS code (based on 
the UPS code) has been developed to compute MHD 
flows in the low magnetic Reynolds number regime. 
The MHD effects are modeled by introducing the ap- 
propriate source terms into the PNS equations. Up- 
stream elliptic effects cam be accounted for by us- 
ing multiple streamwise sweeps with either the IPNS, 
TIPNS, or FBIPNS adgorithms. Turbulence has been 
included by modifying the Baldwin-Lomax turbu- 
lence model [20] to account for MHD effects using 
the approach of Lykoudis [21]. The new code has 
been tested by computing both laminar and turbu- 
lent, supersonic MHD flows over a flat plate. Com- 
parisons have been made with the previous complete 
N-S computations of Dietiker and Hoffmann [18]. In 
addition, the new code has been used to compute the 
supersonic viscous flow inside a rectamgular channel 
designed for MHD experiments [22]. 


( 3 ) 


Ohm's law 

J = (r c (E + VxB) (4) 

where V is the velocity vector, B is the magnetic 
field vector, E is the electric field vector, and J is the 
conduction current density. The flow is assumed to 
be either in chemical equilibrium or in a frozen state. 
The curve fits of Srinivasan et al. [23, 24] are used 
for the thermodynamic and transport properties of 
equilibrium air. 

The governing equations are nondimensionalized 
using the following reference variables. 
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where the superscript * refers to the nondimensional 
quantities. In subsequent sections, the asterisks are 
dropped. 

If the flow variables are assumed to vary in only two 
dimensions (x, y) while the velocity, magnetic, and 
electric fields have components in three dimensions 
(x, y, z), the governing equations can be written in 
the following flux vector form: 


Governing Equations 

The governing equations for a viscous MHD flow with 
a small magnetic Reynolds number are given by [18]: 
Continuity equation 

^ + V - (pV) = 0 (1) 


3U dEi dFi dE v dF„ 
~dt + ~d^ + ~dy = ~dx + ~dy + Smhd 


( 6 ) 


where U is the vector of dependent variables and E, 
and F , are the in viscid flux vectors, and E v and F v 
are the viscous flux vectors. The source term Smhd 
contains all of the MHD effects. The flux vectors are 
given by 


Momentum equation 

d(pV) 


dt 

Energy equation 

djpet) 

dt 


+ V-Lw + p!]=V-f + JxB (2) 
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( 11 ) 


P e t — x/> ( u2 + v 2 + u;2 ) + c — (12) 

2 . 7 - X 

and 7 can be determined from the curve fits of Srini- 
vasan et al. [23] for an equilibrium air flow or is equal 
to a constant (7) for a frozen or perfect gas flow. The 
nondimensional shear stresses and heat fluxes are de- 
fined in the usual manner [11]. 

The governing equations are transformed into com- 
putational space and written in a generalized coordi- 
nate system (£,77) as 

yU t + E f + F, = ^S (13) 

where 

E = (y)(B*-E.)+^(F i -F.) 

F = (y)(B*-B-)+(^)( F<-F w ) (14) 

and J is the Jacobian of the transformation. 

The governing equations are parabolized by drop- 
ping the time derivative term and the streamwise 
direction (£) viscous flow terms in the flux vectors. 
Equation (13) can then be rewritten as 


+ F„ — 


Smhd 


(15) 


where 

E 

F 


7 > + (^' f ‘ 


= (y) (E ,• - K) + (&) (F, - F' v ) (16) 


The prime in the preceding equation indicates that 
the streamwise viscous flow terms have been dropped. 

For turbulent flows, the two-layer Baldwin- Lorn ax 
turbulence model [20] has been modified to account 
for MHD effects. Only the expression for turbulent 
viscosity in the inner layer is changed. This modifi- 
cation for MHD flows is due to Lykoudis [9,21]. 


Numerical Method 


The governing PNS equations with MHD source 
terms have been incorporated into NASA’s upwind 
PNS (UPS) code [13]. These equations can be solved 
very efficiently using a single sweep of the flowfield 
for many applications. For cases where upstream 
(elliptic) effects are important, the flowfield can be 
computed using multiple streamwise sweeps with ei- 
ther the IPNS [14], TIPNS [15], or FBIPNS [16] algo- 
rithms. This iterative process is continued until the 
solution is converged. 

For the iterative PNS (IPNS) method, the E vec- 
tor is split using the Vigneron parameter (u;) [25]. 
This parameter does not need to be changed for the 
present low magnetic Reynolds number formulation. 
In the previous high magnetic Reynolds number code 
[12] it was necessary to modify the Vigneron param- 
eter to account for MHD effects. After splitting, the 
E vector can be written as: 
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(18) 


The streamwise derivative of E is then differenced 
using a forward difference for the “elliptic” portion 

(E p ): 


(f ) j+1 = ift( E «- E . r ) + ( E ?«- E W] 

(19) 
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where the subscript (i 4* 1) denotes the spatial index 
(in the f direction) where the solution is currently 
being computed. The vectors E* +1 and Ef +1 are then 
linearized in the following manner: 

E 'hi = E * + (w),. (u *' +1 ~ Ui) 

E? +1 = E ? + (w) i (U<+ 1 - U<) (20) 


The Jacobians can be represented by 


A' 

A p 


dE‘ 

au 

SEP 

au 


( 21 ) 


After substituting the above linearizations into 
Eq. (19), the expression for the streamwise gradient 
of E becomes 



A* 


(AT-Af)( U, +1 -U,) 


+ (E?+2 — Ef) ] (22) 


The final discretized form of the fluid flow equa- 
tions with MHD source terms is obtained by substi- 
tuting Eq. (22) into Eq. (15) along with the linearized 
expression for the flux in the cross flow plane. The 
final expression becomes: 


Test Case 1: Supersonic laminar and 
turbulent flows over a flat plate with 
applied magnetic field 

In this test case, the supersonic, laminar and turbu- 
lent flow over a flat plate with an applied magnetic 
field is computed. This case corresponds to the flat 
plate case computed previously by Dietiker and Hoff- 
mann [18] using the full N-S equations. A strong 
magnetic field is applied normal to the flow as shown 
in Fig. 1. The dimensional flow parameters for this 
test case are: 


Moo 

= 

2.0 

Poo 

= 

1.076 x 10 s N/m : 

Too 

= 

300 K 

Rcoo 

= 

3.75 x 10 6 

7 

= 

1.4 

L 

= 

0.08 m 

<r e 


800 mho/m 


The plate is assumed to be an adiabatic wall and a 
perfect gas flow is assumed. The magnetic Reynolds 
number (based on the length of the plate) is 0.056 and 
cam be considered negligible when compared to one. 
The normal magnetic field component (B y ) ranges 
in value from 0.0 to 1.2 T. The magnitude of the 
magnetic field cam be represented by the parameter 
m which is defined [18] by 



= RHS (23) 

where 

RHS = _J_[ (E f + 2 )*-(E?) fe+1 ] 

and the superscript k-fl denotes the current iteration 
(i.e. sweep) level. 


Numerical Results 

In order to investigate the utility and accuracy of 
the present PNS approach of solving MHD flowfields 
at low magnetic Reynolds numbers, a few basic test 
cases were computed. The supersonic viscous flow in 
these cases is altered by the presence of the magnetic 
and electric fields which are applied to the flow. 


Poo Uoo 


(24) 


and has units of (1/m). For = 1.2 T, m is equal 
to 1.33. 

A highly stretched grid consisting of 50 points in 
the normal direction was used to compute this case. 
The first point off the wall was located at 2 x 10” 7 
m. Initially, the flow was assumed laminar and sev- 
eral values of B y ranging from 0.0 (no magnetic field) 
to 1.2 T were used. The velocity and temperature 
profiles at x = 0.06 m are shown in Figs. 2 and 3 for 
B y = 0.0 T, 1.0 T, and 1.2 T. The velocity profiles 
are compared to the N-S results of Dietiker and Hoff- 
mann in Fig. 2 and show excellent agreement. The 
magnetic field generates a Lorentz force which acts in 
a direction opposite to the flow. Thus, the flow is de- 
celerated as the magnetic field is increased as seen in 
Fig. 2. For B y — 1.2 T the flow is slightly separated. 
The temperature profiles cannot be compared at this 
time since no temperature data is given in Ref. [18]. 

The turbulent flow over the flat plate was then 
computed using the modified Baldwin-Lomax turbu- 
lence model that accounts for MHD effects. The flow 
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was assumed laminar prior to the point (z = 0.04 m) 
where transition from laminar to turbulent flow was 
triggered. Again, several values of B y ranging from 
0.0 to 1.2 T were used in the computations. The tur- 
bulent velocity and temperature profiles at z = 0.06 
m are shown in Figs. 4 and 5 for B y = 0.0, 1.0 T, 
and 1.2 T. The turbulent velocity profiles in Fig. 4 
are in good agreement with the results of Ref. [18]. 
The variation of skin friction coefficient is shown in 
Fig. 6. The present laminar/ turbulent skin friction 
variations are compared with the results of Ref. [18] 
and show' good agreement. The difference in results 
near the transition point may be due to the coarse 
grid and smoothing used in Ref. [18]. 

All of the present laminar computations were per- 
formed using a single sweep of the flowfield except for 
the separated flow case ( B y — 1.2 T). For this case 
as well as for all the turbulent cases, multiple sweeps 
were used to account for upstream effects. 

Test Case 2: Supersonic viscous flow in 
a rectangular MHD accelerator 

In this test case, the supersonic flow in an experi- 
ment ad MHD channel is simulated. This facility is 
currently being built at NASA Ames Research Center 
by D. W. Bogdanoff, C. Park, and U. B. Mehta [22] 
to study criticad technologies related to MHD bypass 
scramjet engines. The channel is about a half meter 
long and contains a nozzle section, a center section, 
and an aiccelerator section. The channel has a uni- 
form width of 2.03 cm. Magnetic and electric fields 
can be imposed upon the flow in the aiccelerator sec- 
tion. A schematic of the MHD accelerator section is 
shown in Fig. 7. 

This test case was previously computed by 
R. W. MacCormack [10] using the full N-S equations 
coupled with the electromagnetodynamic equations. 
The electrical conductivity in his calculations was set 
at 1.0 x 10 5 mho/m resulting in a very large magnetic 
Reynolds number. In the present study, the calcu- 
lations are performed in the low magnetic Reynolds 
number regime using a realistic value of electrical con- 
ductivity. The flow is computed in two dimensions, 
but later will be extended to three dimensions. Be- 
cause of flow symmetry, only half of the channel is 
computed in the 2-D calculations. 

The flow in the nozzle section and the center sec- 
tion was computed using a combination of the OVER- 
FLOW code [26] and the present PNS code (without 
MHD effects). The initial conditions for the nozzle 
(flow at rest) were: 

po = 8.0 x 10 5 N/m 2 


T 0 = 7500 K 

The laminar flow was assumed to be in chemical equi- 
librium. The computed flowfield at the end of the 
center section was then used as the starting solu- 
tion for the flow calculation of the accelerator section. 
The MHD parameters used in the accelerator section 
were: 



= 50 mho/m 

By 

= 1.5 T 

E z 

= —Ku e B y 

i?e m 

= 0.05 


where the load factor (if) ranged in values from 0.0 to 
1.4, and the centerline velocity (u c ) at the beginning 
of the accelerator section had a value of 3162 m/s. 

The velocity profiles at the end of the accelerator 
section sure shown in Fig. 8 for different load factors. 
The velocity profile with no electric or magnetic fields 
is denoted by K = 0. The increase in the centerline 
velocity with distance (z) for various load factors is 
shown in Fig. 9. The centerline velocity increases 
by about 30% with a load factor of 1.4. It should 
be noted that the flow decelerates because of friction 
when no electric or magnetic fields are applied. 

Concluding Remarks 

In this study, a new parabolized Navier-Stokes al- 
gorithm has been developed to efficiently compute 
MHD flows in the low magnetic Reynolds number 
regime. The new algorithm has been used to com- 
pute both laminar and turbulent, supersonic, MHD 
flows over flat plates and in a rectangular accelera- 
tor section. Although only limited results have been 
obtained thus far, it can be seen that the present 
approach is quite promising. Computations of other 
test cases are currently underway in order to validate 
the current method. 

Acknowledgments 

This work was supported by NASA Ames Research 
Center under Grants NCC2-5379 and NCC2-5517 
and by Iowa State University. The Technical Mon- 
itor for the NASA grant is Dr. Unmeel B. Mehta. 

References 

[1] Gaitonde, D. V., “Development of a Solver for 3- 
D Non-Ideal Magnetogasdynamics AIAA Pa- 
per 99-3610, June 1999. 


5 

American Institute of Aeronautics and Astronautics 


[2] Damevin, H. M., Dietiker, J.-F., and Hoffmann, 
K. A., “Hypersonic Flow Computations with 
Magnetic Field,’’ AIAA Paper 2000-0451, Jan. 
2000 . 

[3] Hoffmann, K. A., Damevin, H. M., and Di- 
etiker, J.-F , “Numerical Simulation of Hyper- 
sonic Magnetohydrodynamic Flows,” AIAA Pa- 
per 2000-2259, June 2000. 

[4] Gaitonde, D. V. and Poggie, J., “Simulation of 
MHD Flow Control Techniques,” AIAA Paper 
2000-2326, June 2000. 

[5] MacCormack, R. W., “Numerical Computa- 
tion in Magnetofluid Dynamics,” Computational 
Fluid Dynamics for the 21st Century , Kyoto, 
Japan, July 2000. 

[6] Deb, P. and Agarwal, R., “Numerical Study of 
MHD-Bypass Scramjet Inlets with Finite-Rate 
Chemistry,” AIAA Paper 2001-0794, Jan. 2001. 

[7] MacCormack, R. W., “A Computational 
Method for Magneto-Fluid Dynamics,” AIAA 
Paper 2001-2735, June 2001. 

[8] Gaitonde, D. V. and Poggie, J., “An Implicit 
Technique for 3-D Turbulent MGD with the 
Generalized Ohm’s Law,” AIAA Paper 2001- 
2736, June 2001. 

[9] Dietiker, J.-F. and Hoffmann, K. A., “Numerical 
Simulation of Turbulent Magnetohydrodynamic 
Flows,” AIAA Paper 2001-2737, June 2001. 

[10] MacCormack, R. W., “Three Dimensional 
Magneto-Fluid Dynamics Algorithm Develop- 
ment,” AIAA Paper 2002-0197, Jan. 2002. 

[11] Tannehill, J. C., Anderson, D. A., and Pletcher, 

R. H., Computational Fluid Mechanics and Heat 
Transfer , Taylor and Francis, Washington, D.C., 
1997. 

[12] Kato, H., Tannehill, J. C., Ramesh, M. D., and 
Mehta, U. B., “Computation of Magnetohydro- 
dynamic Flows Using an Iterative PNS Algo- 
rithm,” AIAA Paper 2002-0202, Jan. 2002. 

[13] Lawrence, S. L., Tannehill, J. C., and Chaussee, 
D. S., “Upwind Algorithm for the Parabo- 
lized Navier-Stokes Equations,” AIAA Journal , 
Vol. 27, No. 9, Sept. 1989, pp. 1975-1983. 

[14] Miller, J. H., Tannehill, J. C., and Lawrence, 

S. L., “PNS Algorithm for Solving Supersonic 
Flows with Upstream Influences,” AIAA Paper 
98-0226, Jan. 1998. 


[15] Tannehill, J. C., Miller, J. H., and Lawrence, 

S. L., “Development of an Iterative PNS Code 
for Separated Flows,” AIAA Paper 99-3361, 
June 1999. 

[16] Kato, H. and Tannehill, J. C., “Development 
of a Forward-Backward Sweeping Parabolized 
Navier-Stokes Algorithm,” AIAA Paper 2002- 
0735, Jan. 2002. 

[17] Munipalli, R., Anderson, D. A., and Kim, H., 
“Two-Temperature Computations of Ionizing 
Air with MHD Effects,” AIAA Paper 2000-0450, 
Jan. 2000. 

[18] Dietiker, J.-F. and Hoffmann, K. A., “Boundary 
Layer Control in Magnetohydrodynamic Flows,” 
AIAA Paper 2002-0130, Jan. 2002. 

[19] Cheng, F., Zhong, X., Gogineni, S., and Kimmel, 
R. L., “Effect of Applied Magnetic Field on the 
Instability of Mach 4.5 Boundary Layer over a 
Flat Plate,” AIAA Paper 2002-0351, Jan. 2002. 

[20] Baldwin, B. S. and Lomax, H., “Thin Layer 
Approximation and Algebraic Model for Sep- 
arated Turbulent Flows,” AIAA- Paper 78-257, 
Jan. 1978. 

[21] Lykoudis, P. S., “Magneto Fluid Mechanics 
Channel Flow, II Theory,” The Physics of Flu- 
ids , Vol. 10, No. 5, May 1967, pp. 1002-1007. 

[22] Bogdanoff, D. W., Park, C., and Mehta, U. B., 
“Experimental Demonstration of Magneto- 
Hydrodynamic (MHD) Acceleration - Facility 
and Conductivity Measurements,” NASA TM- 
2001-210922, July 2001. 

[23] Srinivasan, S., Tannehill, J. C., and Weilmuen- 
ster, K. J., “Simplified Curve Fits for the 
Thermodynamic Properties of Equilibrium Air,” 
NASA RP 1181, Aug. 1987. 

[24] Srinivasan, S., Tannehill, J. C., and Weilmuen- 
ster, K. J., “Simplified Curve Fits for the Trans- 
port Properties of Equilibrium Air,” NASA CR 
178411, 1987. 

[25] Vigneron, Y. C., Rakich, J. V., and Tan- 
nehill, J. C., “Calculation of Supersonic Flow 
over Delta Wings with Sharp Subsonic Leading 
Edges,” AIAA Paper 78-1137, July 1978. 

[26] Buning, P. G., Jespersen, D. C., and Pulliam, 

T. H., OVERFLOW Manual , Version 1.7v , 
NASA Ames Research Center, Moffett Field, 
California, June 1997. 


6 

American Institute of Aeronautics and Astronautics 









8 

American Institute of Aeronautics and Astronautics 







x, m 


Figure 6: Laminar/ turbulent skin friction coefficient 
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Figure 9: Centerline velocity for various load factors 
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Abstract 

The 3-D supersonic viscous flow in an experimental 
MHD channel has been numerically simulated. The 
experimental MHD channel is currently in operation 
at NASA Ames Research Center. The channel con- 
tains a nozzle section, a center section, and an ac- 
celerator section where magnetic and electric fields 
can be imposed on the flow. In recent tests, veloc- 
ity increases of up to 40% have been achieved in the 
accelerator section. The flow in the channel is numer- 
ically computed using a new 3-D parabolised Navier- 
Stokes (PNS) algorithm that has been developed to 
efficiently compute MHD flows in the low magnetic 
Reynolds number regime. The MHD effects are mod- 
eled by introducing source terms into the PNS equa- 
tions which can then be solved in a very efficient man- 
ner. To account for upstream (elliptic) effects, the 
flowfield can be computed using multiple streamwise 
sweeps with an iterated PNS algorithm. The new al- 
gorithm has been used to compute two test cases that 
match the experimental conditions. In both cases, 
magnetic and electric fields are applied to the flow. 
The computed results are in good agreement with the 
available experimental data. 

Introduction 

Magnetohyrodynamics (MHD) can be utilized to im- 
prove performance and extend the operational range 
of many systems. Potential applications include hy- 
personic cruise, advanced Earth-to-orbit propulsion, 
chemical and nuclear space propulsion, regenerative 
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aerobraking, onboard flow control systems, test fa- 
cilities, launch assist, and power generation. One of 
the critical technologies associated with these appli- 
cations is MHD acceleration. In order to study MHD 
acceleration, an experimental MHD channel has been 
built at NASA Ames Research Center by D. W. Bog- 
danoff, C. Park, and U. B. Mehta [1,2]. The channel 
is about a half meter long and contains a nozzle sec- 
tion, a caiter section, and an accelerator section. The 
channel has a uniform width of 2.03 cm. Magnetic 
and electric fields can be imposed upon the flow in 
the accelerator section. A cross section of the MHD 
ch ann el is shown in Fig. 1. 

In the present study, the flow in the experimen- 
tal MHD channel is numerically simulated. Flow- 
fields involving MHD effects have typically been 
computed [3-15] by solving the complete Navier- 
Stokes (N-S) equations for fluid flow in conjunction 
with Maxwell’s equations of electromagnetodynam- 
ics. When chemistry and turbulence effects are also 
included, the computational effort required to solve 
the resulting coupled system of partial differential 
equations is extremely formidable. One possible rem- 
edy to this problem is to use the parabolized Navier- 
Stokes (PNS) equations in place of the N-S equar- 
tions. The PNS equations can be used to compute 
three-dimensional, supersonic viscous flowfields in a 
very efficient manner [16]. This efficiency is achieved 
because the equations can be solved using a space- 
marching technique as opposed to the time-marching 
technique that is normally employed for the complete 
N-S equations. 

Recently, the present authors have developed PNS 
codes to solve 2-D supersonic MHD flowfields in both 
the high and low magnetic Reynolds number regimes 
[17,18]. The magnetic Reynolds number is defined 
as Rc-rn — <7 e p c Foo£ where a c is the electrical con- 
ductivity, fj> c is the magnetic permeability, Fqo is the 
freestream velocity, and L is the reference length. 
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The new MUD PNS codes are based on NASA’s up- 
wind PNS (UPS) code which was originally devel- 
oped by Lawrence et al. [19]- The UPS code solves 
the PNS equations using a fully conservative, finite- 
volume approach in a general nonorthogonal coordi- 
nate system. The UPS code has previously been ex- 
tended to permit the computation of flowfields with 
strong upstream influences- In regions where strong 
upstream influences are present, the governing equa- 
tions are solved using multiple sweeps (i.e. itera- 
tions). As a result of this approach, a complete flow- 
field can be computed more efficiently (in terms of 
computer time and storage) than with a standard N- 
S solver which marches the entire solution in time. 
Three iterative PNS algorithms called IPNS, TIPNS, 
and FBIPNS have been developed and are described 
in Refe. [20-22]. 

For many aerospace applications, including the 
present experimental MHD channel, the electrical 
conductivity of the fluid is low and hence the mag- 
netic Reynolds number is small. In these cases, it 
makes sense to use the low magnetic Reynolds num- 
ber assumption and reduce the complexity of the gov- 
erning equations. The MHD effects are modeled with 
the introduction of source terms into the fluid flow 
equations* as was done in the present low magnetic 
Reynolds number PNS code [18]. This code has been 
extended to three-dimensions in the present study 
and is used for the numerical simulations of the flow 
in the experimental MHD channel. 


Governing Equations 

The governing equations for a viscous MHD flow 
with a small magnetic Reynolds number are given 
by [14]: 

Continuity equation 

|? + V(pV) = 0 (1) 

Momentum equation 

^+V-[pW+pi]=V-f + JxB (2) 
Energy equation 

^p+V-[(pe* +1>) V] = V-(V-f)— V-U+B-J (3) 
Ohm’s law 

J = cr e (E + VxB) (4) 

where V is the velocity vector, B is the magnetic 
field vector, E is the electric field vector, and J is the 
conduction current density. 


The governing equations are nondimensi onalized 
using the following reference variables: 




L ’ 


U*jV*,W* 


u, v, w 


t* 


Uoot 


L 


P = 
= 


P 

Poo 

JL 

Ul 


T* = 


P = 


fL 

Poo u* 


P = 


PooUl 
P 


Poo 


(5) 


R* R* R* 1 ^ ^ ip* zp* ip* 

Uooy/Pepoo ^oo^/pepo o 


Pc = 


Pc 

Pc oo 


= 1, *; = 


where the superscript * refers to the nondimen- 
sional quantities. For convenience, the asterisks are 
dropped in the following equations. 

The governing equations written in vector form in 
a 3-D Cartesian coordinate system become 


OTJdEidTidGi _ dE v dF„dG v i c 
dt + dx + dy + dz dx + dy + dz +Smhd 

( 6 ) 

where U is the vector of dependent variables, E*, F » 
and G{ are the in viscid flux vectors, and E w , F„ and 
G v are the viscous flux vectors. The source term 
Smhd contains all of the MHD effects. The flux vec- 
tors are given by 


U = [ p, pu, pv, pw, pet ] T (7) 


E t 


pu 


pv 

pu 2 +p 


puv 

puv 

, Fi = 

pro 1 + p ' 

puw 


pvw 

. (pe t +p)u . 


. (pe t -hp)u _ 


pw 


G i = 


puw 
pv w 
pw 2 +p 

(pe t +p)w 


( 8 ) 


E„ 


F„ 


t IT XX + VT xy + WT XZ - q x J 


* yx 

T yy 
r yz 

UTyx + VTyy + XVTy Z — J 


(9) 


( 10 ) 
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(11) 


G„ = 


0 

7 ’** 

T ZV 

Tzz 

UTzx + VTzy + WT zz - 


SMHD = -RCm 


0 

B Z (E V + wB x — uB z ) 
—B V (E Z 4 - uB y — vB x ) 
B x {E z + uj By — vB x ) 
-B Z (E X + vBz — wB y ) 
By(E x + vB x — tvBy) 
-BxjE^ + wB x — uB z ) 
E X (E X + vB z — wBy) 
+Ey(Ey + wB x — u B z ) 
+E Z {E Z + u By — vB z ) 


where 


pet = ip^ + t^ + ur 2 ) + 


(12) 


(13) 


and the nondimensional shear stresses and heat fluxes 
are defined in the usual manner [16]- 

The flow can be computed assuming either a con- 
stant 7 or by using the simplified curve fits of Srini- 
vasan et al. [23,24] for the thermodynamic and trans- 
port properties of equilibrium air. For the latter case, 
7 in Equation (13) is determined using the curve fit 
7 = 7(e,p). In future calculations, the flow will be 
computed in chemical nonequilibrium. 

The governing equations are transformed into com- 
putational space and written in a generalized coordi- 
nate system (f , rj, C) as 


iu t + E e +F„ + G c = ^p (14) 

where 

E = (y) (E< - E„) + (&) (f 4 - F.) 

+ ( 7 ) (g ‘- g ’> 

F = (jj) (E» — E«) + (jj) (Fi — F„) 

+ { T j) (G<_G " ) 

G = ^(E,-E.)+0^(F 4 -F.) 

+ ^y^ ( G « ~ Gw ) (15) 

and J is the Jacobian of the transformation. 

The governing equations are parabolized by drop- 
ping the time derivative term and the streamwise 


direction (f) viscous flow terms in the flux vectors. 
Equation (14) can then be rewritten as 

E ? + F„ + G c = ^p (16) 

where 

* - (tt) ^ + (!f) Fi + ( /) G< 

F = (^)(B 1 -E'J + (i)(F < -r v ) 

+ (t) (G<_g;) 

g = ( Ei -E;) + ^) (Fi-r v ) 

+ (y) (G < - G t) (17) 

The primes in the preceding equations indicate that 
the streamwise viscous flow terms have been dropped. 

For turbulent flows, the two-layer Baldwin-Lomax 
turbulence model [25] has been modified to account 
for MHD effects. Only the expression for turbulent 
viscosity in the inner layer is changed. This modifi- 
cation for MHD flows is due to Lykoudis [26]. 


Numerical Method 

The governing PNS equations with MHD source 
terms have been incorporated into NASA’s upwind 
PNS (UPS) code [19]. These equations can be solved 
very efficiently using a single sweep of the flowfield 
for many applications. For cases where upstream 
(elliptic) effects are important, the flowfield can be 
computed using multiple streamwise sweeps with ei- 
ther the IPNS [20], TIPNS [21], or FBIPNS [22] algo- 
rithms. This iterative process is continued until the 
solution is converged. 

For the iterative PNS (IPNS) method, the E vec- 
tor is split using the Vigneron parameter (uj) [27]. 
This parameter does not need to be changed for the 
present low magnetic Reynolds number formulation. 
In the previous high magnetic Reynolds number code 
[17] it was necessary to modify the Vigneron param- 
eter to account for MHD effects. After splitting, the 
E vector can be written as: 

E = E*+E p (18) 

where 


pu 


pv 

pu 2 + up 


puv 

puv 

+ i 

pv 2 +ujp 

puw 


pvw 

_ (pe t +p)u _ 


_ (pe t +p)v _ 
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pw 

puw 

pvw 

pCJ 2 + up 
(pe t +p)w 


p 

0 

(l-u)p 

c 

0 

0 

w = — 

J 

0 

0 

0 

+ 7 

(1 -u)p 
0 
0 



0 

0 

0 

(l-w)p 

0 


(19) 


The streamwise derivative of E is then differenced 
using a backward difference for E* and a forward dif- 
ference for the “elliptic” portion (E J> ): 


(^) = ^ [(®<+i - E*) + (Ef +2 - Ef +1 )] 

( 20 ) 

where the subscript (t 4- 1) denotes the spatial index 
(in the f direction) where the solution is currently 
being computed. The vectors E* +1 and E^ are then 
linearized in the following manner: 


E; +1 = E*+^^(U, +1 -U t ) 

Ef+i = Ef+(^) j (U i +i-U t ) 
The Jacobians can be represented by 


( 21 ) 


A* 

A p 


3E* 

au 

dW 
d U 


( 22 ) 


After substituting the above linearizations into 
Eq. (20), the expression for the streamwise gradient 
of E becomes 


+ (Ef+2 ~ Ef) ] (23) 


The final discretized form of the fluid flow equa- 
tions with MHD source terms is obtained by substi- 
tuting Eq. (23) into Eq. (16) along with the linearized 
expressions for the fluxes in the cross flow plane. The 
final expression becomes: 


±(a--a?)+ 9 < dF 


.d_(9G 

• r\ > 


-I fc+1 


dv \dv j, ac 

= RHS 


(AUi)* +1 

(24) 


where 


(AUi)* +1 = (U i+1 -U i )‘ +1 



and the superscript k+1 denotes the current itera- 
tion (i.e. sweep) level. In the preceding equation, the 
MHD source term, Smhd, is treated explicitly since it 
is evaluated using the velocity at station i (V*). For 
most cases, this will not degrade the accuracy of the 
solution since Af is small and the velocity changes 
slowly. If this is not the case, a predictor-corrector 
procedure can be implemented whereby a predicted 
velocity at station i+1 (V\~ ) is first obtained using 
Eq. (24). The solution at station i-Hl is then recom- 
puted by evaluating Smhd with V~ - 


Numerical Results 

The numerical calculation of the supersonic flow in 
the experimental MHD channel is now discussed. 
This flowfield was previously computed by MacCor- 
mack [12] using the full N-S equations coupled with 
the electromagnetodynamic equations. The electrical 
conductivity in his calculations was set at 1.0 x 10 4 5 
mho/m resulting in a very large magnetic Reynolds 
number. In the present study, the calculations are 
performed in the low magnetic Reynolds number 
regime using a realistic value of electrical conductiv- 
ity. Both 2-D and 3-D results have been obtained. 

The flow in the nozzle section was computed using 
a combination of the OVERFLOW code [28] and the 
present PNS code (without MHD effects). For the 3- 
D OVERFLOW nozzle calculation, a highly stretched 
grid consisting of 130 x 50 x 50 grid points was used. 
The normal grid spacing at the wall was 1.0 x 10“ 5 
m. For the PNS calculation of the flow in the re- 
mainder of the nozzle and the rest of the MHD chan- 
nel, a highly stretched grid consisting of 60 points in 
both the y and z directions was used and the nor- 
mal grid spacing at the wall was 1.0 x 10“ 5 m. As 
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a consequence of flow symmetry, only one-fourth of 
the channel cross section was computed in the 3-D 
calculations. 

The calculations were performed assuming turbu- 
lent flow throughout the MHD channel. In addition 
the flow was computed using a constant 7 of 1.25 
to simulate equilibrium air. The channel wall tem- 
perature was assumed to be isothermal since steady 
flow conditions were maintained in the experiment for 
only about 1.2 milliseconds. A schematic of the pow- 
ered portion of the MHD channel with the directions 
of the applied magnetic and electric fields is shown 
in Fig. 2. The values of the electrical conductivity 
(cr e ), the magnetic field (jB z ), and the electric field 
(E y ) were kept constant in the powered portion of 
the channel. Two test cases corresponding to Runs 
15 and 16 of the NASA Ames experiments [29] were 
computed in this study and are now discussed. 

Test Case 1: NASA Ames MHD Run 15 

(Vcap. = 32QV) 

The dimensional flow parameters for this test case 
are: 

Po = 9.10 x IQ 5 N/m 2 

T 0 = 5560 K 

T w = 300 K 

<r e ~ 130 mho/m 

B z = 0.0, 0.92 T 

E y = 0, 3955, 5000, 6000 Y/m 

where the subscript o denotes total conditions at the 
nozzle entrance and w denotes wall conditions. 

This case was computed using several different elec- 
tric field strengths in order to properly simulate the 
experiment. In the experiment, the voltage applied 
to the electrodes was approximately 134 V for this 
case, however, due to the sheath voltage drop, the 
actual voltage applied to the flow is smaller than the 
electrode voltage. The voltage drop was measured for 
the central inviscid core flow, and was approximately 
67 V [2]. Since the boundary layer is computed in the 
numerical solution, the applied electric field must be 
approximately the voltage drop across the electrodes 
minus the sheath voltage drop. Unfortunately, it is 
not a trivial task to measure the sheath voltage drop. 
Therefore, several different electric fields were cho- 
sen in the numerical calculations so that the corre- 
sponding voltage drop across the electrodes would be 
between 67 V and 134 V. The voltage drop of 67 V 
corresponds to E y = 3955 V/m at the center of the 
accelerator section and a voltage drop of 101.6 V cor- 
responds to E y = 6000 V/m, with a sheath voltage 
drop of 32.4 V. 


The computed streamwise variation of static pres- 
sure for the 2-D calculations is shown in Fig. 3 for the 
different electric field strengths. The pressure varia- 
tion with no electric field or magnetic field is denoted 
by E y = 0. The results are in reasonable agreement 
with the experiment. The numerical results show an 
increase in static pressure as the electric field strength 
is increased. 

The computed streamwise variation of averaged ve- 
locity for the 2-D calculations is shown in Fig. 4. The 
velocities are averaged across the channel cross sec- 
tion and normalized using the entrance velocity to be 
consistent with the experiment. In the experiment, 
the velocities were obtained by measuring the volt- 
age generated by the flow at the last electrode pair 
(19) which is unpowered. This procedure inherently 
involves an averaging of the velocity profile. The nu- 
merical results indicate an increase in the averaged 
velocity of about 26% with E y = 6000 V/m and this 
compares well with the experimental value of approx- 
imately 27%. The velocity vector plots for Ey = 0 
and E y — 6000 V /m are shown in Fig. 5. 

The computed streamwise variation of averaged ve- 
locity for the 3-D calculations is shown in Fig. 6. The 
velocities are averaged and normalized in the same 
manner as for the 2-D calculations- The numerical 
results indicate an increase in the averaged velocity 
of about 26% with E y = 6000 V/m which is the same 
value obtained with the 2-D calculations. 

Test Case 2: NASA Ames MHD Run 16 

(V cap . = 380V) 

The dimensional flow parameters for this test case 
are: 

Po = 9.92 x10 s N/m 2 
To = 5560 K 
T w — 300 K 
a e = 140 mho/m 
B z = 0.0, 0.92 T 

E y = 0, 4309, 5000, 6000, 7000 V/m 

This test case was also computed using several dif- 
ferent electric field strengths in order to properly 
simulate the experiment. The computed stream- 
wise variation of static pressure for the 2-D calcula- 
tions is shown in Fig. 7 for the different electric field 
strengths. The computed pressures are in good agree- 
ment with the experimental pressures. The computed 
streamwise variation of averaged velocity for the 2-D 
calculations is shown in Fig. 8. The numerical results 
indicate an increase in the averaged velocity of about 
31% with E y = 7000 V /m. This is less than the value 
of 38% that was obtained in the experiment. The cor- 
responding streamwise variation of averaged velocity 
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for the 3-D calculations is shown in Fig. 9. The nu- 
merical results indicate an increase in the averaged 
velocity of about 30% with E y = 7000 V /m which is 
again less than the experimental value of 38%. 

Concluding Remarks 

In this study, a new 3-D parabolized Navier-Stokes 
algorithm has been developed to efficiently compute 
MHD flows in the low magnetic Reynolds number 
regime. The new algorithm has been used to com- 
pute the flow in the NASA Ames experimental MHD 
channel for Runs 15 and 16. The numerical results 
are in good agreement with most of the experimental 
results. 
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Figure 2: Schematic of powered portion of MHD channel 
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Figure 4: Streamwise variation of averaged velocity (2-D) 
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Figure 6: Streamwise variation of averaged velocity (3-D) 
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Figure 7: Streamwise variation of static pressure (2-D) 



Figure 8: Streamwise variation of averaged velocity (2-D) 
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Figure 9: Stream wise variation of averaged velocity (3-D) 
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Abstract 

The 3-D nonequilibrium seeded air flow in the NAS A- 
Ames experimental MHD channel has been numeri- 
cally simulated. The channel contains a nozzle sec- 
tion, a center section, and an accelerator section 
where magnetic and electric fields can be imposed 
on the flow. In recent tests, velocity increases of up 
to 40% have been achieved in the accelerator section. 
The flow in the channel is numerically computed us- 
ing a 3-D parabolized Navier-Stokes (PNS) algorithm 
that has been developed to efficiently compute MHD 
flows in the low magnetic Reynolds number regime. 
The MHD effects are modeled by introducing source 
terms into the PNS equations which can then be 
solved in a very efficient manner. The algorithm has 
been extended in the present study to account for 
nonequilibrium seeded air flows. The electrical con- 
ductivity of the flow is determined using the program 
of Park. The new algorithm has been used to com- 
pute two test cases that match the experimental con- 
ditions. In both cases, magnetic and electric fields 
are applied to the seeded flow. The computed results 
are in good agreement with the experimental data. 

Introduction 

Magnetohyrodynamics (MHD) can be utilized to im- 
prove performance and extend the operational range 
of many systems. Potential applications include hy- 
personic cruise, advanced Earth-to-orbit propulsion, 
chemical and nuclear space propulsion, regenerative 
aerobraking, onboard flow control systems, test fa- 
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cilities, launch assist, and power generation. One of 
the critical technologies associated with these appli- 
cations is MHD acceleration. In order to study MHD 
acceleration, an experimental MHD channel has been 
built at NASA Ames Research Center by D. W. Bog- 
danoff, C. Park, and U. B. Mehta [1,2]. The channel 
is about a half meter long and contains a nozzle sec- 
tion, a center section, and an accelerator section. The 
channel has a uniform width of 2.03 cm. Magnetic 
and electric fields can be imposed upon the flow in 
the accelerator section. A cross section of the MHD 
channel is shown in Fig. 1. 

In the present study, the flow in the experimen- 
tal MHD channel is numerically simulated. Flow- 
fields involving MHD effects have typically been 
computed [3-15] by solving the complete Navier- 
Stokes (N-S) equations for fluid flow in conjunction 
with Maxwell’s equations of electromagnetodynam- 
ics. When chemistry and turbulence effects are also 
included, the computational effort required to solve 
the resulting coupled system of partial differential 
equations is extremely formidable. One possible rem- 
edy to this problem is to use the parabolized Navier- 
Stokes (PNS) equations in place of the N-S equa- 
tions. The PNS equations can be used to compute 
three-dimensional, supersonic viscous flowfields in a 
very efficient manner [16]. This efficiency is achieved 
because the equations can be solved using a space- 
marching technique as opposed to the time-marching 
technique that is normally employed for the complete 
N-S equations. 

Recently, the present authors have developed PNS 
codes to solve 2-D and 3-D supersonic MHD flowfields 
in both the high and low magnetic Reynolds number 
regimes [17-19]. The magnetic Reynolds number is 
defined as Re m = a e ^VooL where a e is the electrical 
conductivity, /x e is the magnetic permeability, is 
the ffeestream velocity, and L is the reference length. 
The new MHD PNS codes are based on NASA’s up- 
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wind PNS (UPS) code which was originally developed 
by Lawrence et al. [ 20 ]. The UPS code solves the PNS 
equations using Roe’s scheme in a fully conservative, 
finite- volume approach in general nonorthogonal co- 
ordinates. 

For many aerospace applications, including the 
present experimental MHD channel, the electrical 
conductivity of the fluid is low and hence the mag- 
netic Reynolds number is small. In these cases, it 
makes sense to use the low magnetic Reynolds num- 
ber assumption and reduce the complexity of the gov- 
erning equations. The MHD effects are modeled with 
the introduction of source terms into the PNS equar 
tions. 

Previously [19], the present authors used the low 
magnetic Reynolds PNS code to compute both 2 - 
D and 3-D flows in the NASA- Ames MHD channe l. 
These perfect gas (7 = 1.25) calculations assumed 
that the magnetic and electric fields, as well as the 
electrical conductivity, were constant in the accelera- 
tor section. In the present study, the 3-D simulations 
have been extended to include both equilibrium air 
flows as well as nonequilibrium seeded air flows. For 
the latter case, the electrical conductivity is variable 
and is computed using the program of Park [21]. 


Governing Equations 

Magnetogasdynamic Equations 

The governing equations for a viscous MHD flow 
with a small magnetic Reynolds number are given 
by [14]: 

Continuity equation 

§+V-(pV )=0 ( 1 ) 

Momentum equation 

^ + V- [pW + plj =V-f.+ JxB ( 2 ) 


The governing magnetogasdynamic equations are 
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where the superscript * refers to the nondimen- 
sional quantities. For convenience, the asterisks are 
dropped in the following equations. 

The governing equations written in vector form in 
a 3-D Cartesian coordinate system become 


dU dE* ( dF i t dG i dE v dF v dG v B 
~dT + = ~dx + ~dy’ + ~dT +SuHD 

( 6 ) 

where U is the vector of dependent variables, E it F* 
and G» are the inviscid flux vectors, and E v , F v and 
G v are the viscous flux vectors. The source term 
Smhd contains all of the MHD effects. The flux vec- 
tors are given by 
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Energy equation 
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Ohm’s law 

J = <7 e (E + VxB) (4) 

where V is the velocity vector, B is the magnetic 
field vector, E is the electric field vector, and J is the 
conduction current density. 
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where 
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and the nondimensional shear stresses and heat fluxes 
are defined in the usual manner [16]. 

The governing equations are transformed into com- 
putational space and written in a generalized coordi- 
nate system (f,7j,C) as 


jU t + E e + F„ + G c = ^p ( 14 ) 

where 
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and J is the Jacobian of the transformation. 

The governing equations are parabolized by drop- 
ping the time derivative term and the streamwise 
direction (f) viscous flow terms in the flux vectors. 
Equation (14) can then be rewritten as 
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The primes in the preceding equations indicate that 
the streamwise viscous flow terms have been dropped. 

For turbulent flows, the two-layer B aldwin- Lomax 
turbulence model [22] has been modified to account 
for MHD effects. Only the expression for turbulent 
viscosity in the inner layer is changed. This modifi- 
cation for MHD flows is due to Lykoudis [23]. 

In order to “close” the preceding system of PNS 
equations, relations between the thermodynamic 
variables are required along with expressions for the 
transport properties /x and k. For a perfect gas, the 
pressure is computed from the relation 

p = (j-l)pe (18) 

where 7 = 7oo, and the transport properties are com- 
puted using Sutherland’s formulas [16]. For equilib- 
rium air computations, 7 and all other thermody- 
namic and transport properties are obtained from the 
simplified curve fits of Srinivasan et al. [24,25]. For 
nonequilibrium computations, the thermodynamic 
and transport properties are determined using the 
procedures described in the next section. 


Nonequilibrium Flow Equations 


For nonequilibrium flows, the species continuity equa- 
tions must be solved in addition to the magnetogas- 
dynamic equations given previously The magnet o- 
gasdynamic equations remain the same except for the 
additional term in the energy equation, which is due 
to the diffusion of the species. The nondimensional 
species continuity equations, expressed in 2-D trans- 
formed coordinates for a steady flow, are given by 
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where c s is the mass fraction of species s, u s is the 
nondimensional production term, D is the nondimen- 
sional binary diffusion coefficient, and 0z = P ^ D °° . 
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The chemical model used in the present calcula- 
tions is similar to the clean-air model of Blottner et 
al. [26] and Prabhu et al. [27]. It consists of molecular 
oxygen (0 2 ), atomic oxygen (O), molecular nitrogen 
(iV 2 ), atomic nitrogen (N), nitric oxide {NO), ni- 
tric oxide ion (NO+) and electrons (e“). The follow- 
ing reactions are considered between the constituent 
species. 

(1) 0 2 + M x ^20+Mi 

(2) N 2 + M 2 ^ 2N + M 2 

(3) N 2 + N ^ 2N + N 

(4) i\TO + M 3 ^iV + 0-hM3 

(5) JV'O-hO^Oa + iV' 

(6) N 2 + O^NO + N 

(7) iVO ^ NO* 4* e~ (20) 

where Mi, M 2 , M3, are catalytic third bodies. The 
clean-air chemical model has 7 species (n = 7) and 
seven reactions (m = 7). In order to simulate the 
seeded air flow in the MHD channel, the potassium 
seeding reaction has been added to the above chem- 
istry model. This reaction is the ionization of atomic 
potassium (K) and is given by the following equation: 


(8) K+e~ ^K+ + e- + e~ (21) 

Using the law of mass action, the nondimensional 
mass production rate of species s is 


^ = M,fx, - u’ k<3 ) [K f , k (T) f[]jrr, 
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where 7 r is the nondimensional mole-mass ratio of 
the reactants, M a is the molecular weight of species 
s, i/ ka and u k a are the stoichiometric coefficients and 
th is the number of reactants. Further details on 
the reaction rates and the thermodynamic and trans- 
port properties can be found in Ref. [26]. The electri- 
cal conductivity is determined from the species mole 
fractions, along with the temperature, density, and 
pressure of the gas, using the program of Park [21]. 


Numerical Method 

Solution of PNS Equations 

The governing PNS equations with MHD source 
terms have been incorporated into NASA’s upwind 
PNS (UPS) code [20]. These equations can be solved 
very efficiently using a single sweep of the flowfield 


for many applications. For cases where upstream 
(elliptic) effects awe important, the flowfield can be 
computed using multiple streamwise sweeps with ei- 
ther the IPNS [28], TIPNS [29], or FBEPNS [30] algo- 
rithms. This iterative process is continued until the 
solution is converged. 

For the iterative PNS (IPNS) method, the E vec- 
tor is split using the Vigneron parameter (cj) [31]. 
This parameter does not need to be changed for the 
present low magnetic Reynolds number formulation. 
In the previous high magnetic Reynolds number code 
[17] it was necessary to modify the Vigneron param- 
eter to account for MHD effects. After splitting, the 
E vector can be written as: 


E = E* + EP 
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The streamwise derivative of E is then differenced 
using a backward difference for E* and a forward dif- 
ference for the “elliptic” portion (E p ): 


(f|) = [(E;+1 ’ E:) + (E « ' Er " )1 

(25) 

where the subscript {i + 1) denotes the spatial index 
(in the f direction) where the solution is currently 
being computed. The vectors E* +1 and Ef +1 are then 
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linearized in the following manner: 


Solution of Species Continuity Equations 
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The Jacobians can be represented by 
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After substituting the above linearizations into 
Eq. (25), the expression for the streamwise gradient 
of E becomes 


(f) M - 

+ (Ef+2 - Ef) ] (28) 

The final discretized form of the fluid flow equa- 
tions with MHD source terms is obtained by substi- 
tuting Eq. (28) into Eq. (16) along with the linearized 
expressions for the fluxes in the cross flow plane. The 
final expression becomes: 

[>*-«>+ 1 ( i ), ■ 

+ i(ro)r ,(Au ‘ ) ‘ +i=BHs(29) 

where 

(AU i ) k+1 = (U i+ i-U i ) fe+1 



and the superscript k+1 denotes the current itera- 
tion (i.e. sweep) level. In the preceding equation, the 
MHD source term, Smhd, is treated explicitly since it 
is evaluated using the velocity at station i (V*). For 
most cases, this will not degrade the accuracy of the 
solution since is small and the velocity changes 
slowly. If this is not the case, a predictor-corrector 
procedure can be implemented whereby a predicted 
velocity at station i+1 (Vjj^) is first obtained using 
Eq. (29). The solution at station i+1 is then recom- 
puted by evaluating Smhd with 


For chemical nonequilibrium, the species continuity 
equations, Eq. (19), must be solved in addition to 
the magnetogasdynamic equations. The equations 
have been integrated using the loosely-coupled ap- 
proach of Tannehill et ad. [32]. In this approach, 
the species continuity equations and magnetogasdy- 
namic equations are solved separately. The coupling 
between the two sets of equations is then obtained 
in an approximate manner. The species continuity 
equations are modeled using a second-order-accurate, 
upwind-based TVD scheme for the convective terms 
and second-order-accurate central differences for the 
diffusion terms. The assumption of zero net charge of 
the gas is used to eliminate the electron mass conser- 
vation equation. In addition, the species continuity 
equation for the nth species is eliminated by using 
the requirement that the mass fractions must sum 
to unity. The term representing the rate of produc- 
tion of species, uj 3 , is treated as a source term, and is 
lagged to the previous marching level. 

The coupling between the fluids and the chemistry 
is performed in an approximate manner. First, a fluid 
step is taken from marching station i to i+1 assuming 
frozen chemistry. Then the fluid density and velocity 
at i+1 are used in the solution of the species continu- 
ity equations to obtain species mass fractions at i+1. 
Finally, the species mass fractions, molecular weight 
of mixture, fluid density, and internal energy at i+1 
are used to obtain the new temperature, pressure, 
specific enthalpy, and frozen specific heats at the i+1 
marching station. 

The temperature is obtained by performing a 
Newton-Raphson iteration of the following form: 


Tk+ i_ T k g(T*)-ft 

am 


(30) 


where 

g( T) = J2c 3 h 9 (T) 

5=1 

5 '(t) = y> 3 c p , 4 (T) 

5=1 

and k is the index of iteration. The iterations are 
continued until 


| T k+1 - T k j < € 

where e is a small positive quantity. Once the temper- 
ature is determined, the pressure can be computed us- 
ing Dalton’s law of partial pressures. Further details 
of this procedure can be found in Refe. [32] and [33]. 
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Numerical Results 

The numerical calculation of the 3-D supersonic flow 
in the experimental MHD c h annel is now discussed. 
The flow in the nozzle section was computed using a 
combination of the OVERFLOW code [34] and the 
present PNS code (without MHD effects). For the 
OVERFLOW nozzle calculation, a highly stretched 
grid consisting of 150 x 80 x 80 grid points was used. 
The normal grid spacing at the wall was 1.0 x 10“ 5 
m. For the PNS calculation of the flow in the re- 
mainder of the nozzle and the rest of the MHD chan- 
nel, a highly stretched grid consisting of 90 points in 
both the y and z directions was used and the nor- 
mal grid spacing at the wall was 2.0 x 10~ 6 m. As 
a consequence of flow symmetry, only one-fourth of 
the channel cross section was computed in the 3-D 
calculations. 

The calculations were performed assuming turbu- 
lent flow throughout the MHD channel. The chan- 
nel wall temperature was assumed to be isothermal 
since quasi-steady flow conditions were maintained 
in the experiment for only about 1.2 milliseconds. A 
schematic of the powered portion of the MHD chan- 
nel along with the directions of the applied magnetic 
and electric fields is shown in Fig. 2. The values of the 
magnetic field ( B z ), and the electric field (E v ) were 
kept constant in the powered portion of the channel. 

Three different chemistry models were used in this 
study to simulate the flow in the MHD channel. 
These were: (1) perfect gas (7 = 1.25), (2) equi- 
librium air, and (3) nonequilibrium seeded-air chem- 
istry. For the perfect gas and equilibrium air calcu- 
lations, the electrical conductivity (cr e ) was assumed 
constant. For the nonequilibrium seeded-air calcu- 
lations, the electrical conductivity varied throughout 
the flowfield and was determined using the program 
of Park [21]. The seeding (as in the experiment) con- 
sisted of 1% (by mass) of potassium. Two test cases 
corresponding to Runs 15 and 16 of the NASA Ames 
experiments [2,35] were computed in this study and 
are now discussed. 


Test Case Is 

NASA Ames MHD Rim 15 

C Kap . = 320V) 

The dimensional 
are: 

flow parameters for this test case 

Vo = 

9.10 x 10 s N/m 2 

T 0 - 

5560 K 

T w - 

300 K 

<j e = 

130 mho/m (or variable) 

B z = 

0.0, 0.92 T 

Ey = 

0, 3955, 5000 V/m 


where the subscript o denotes total conditions at the 
nozzle entrance and w denotes wall conditions. 

This case was computed using several different elec- 
tric field strengths in order to properly simulate the 
experiment. In the experiment, the voltage applied 
to the electrodes was approximately 134 V for this 
case, however, due to the sheath voltage drop, the 
actual voltage applied to the flow is smaller than the 
electrode voltage. The voltage drop was measured for 
the central inviscid core flow, and was approximately 
67 V [2]. Since the boundary layer is computed in the 
numerical solution, the applied electric field must be 
approximately the voltage drop across the electrodes 
minus the sheath voltage drop. Unfortunately, it is 
not a trivial task to measure the sheath voltage drop. 
Therefore, several different electric field strengths 
were chosen in the numerical calculations so that 
the corresponding voltage drop across the electrodes 
would be between 67 V and 134 V. The voltage drop 
of 67 V corresponds to E v = 3955 V/m and a voltage 
drop of 84.7 V corresponds to E y = 5000 V/m. 

The computed streamwise variation of static pres- 
sure for the nonequilibrium seeded- air calculations 
is shown in Fig. 3 for the different electric field 
strengths. The pressure variation with no electric 
field or magnetic field is denoted by E v = 0. The re- 
sults for E v = 5000 V/m are in excellent agreement 
with the experiment. The numerical results show an 
increase in static pressure as the electric field strength 
is increased. The computed streamwise variation of 
static pressure for the different chemistry models is 
shown in Fig. 4 for E v = 5000 V/m. The nonequi- 
librium seeded- air model gives the closest agreement 
with the experimental pressures. 

For the nonequilibrium seeded-air computations, 
the electrical conductivity was not constant but var- 
ied throughout the flowfield. The average conduc- 
tivity (averaged over the channel cross section) at 
the center of the powered portion of the channel 
(electrode pair 10) was found to be 130 mho/m for 
E y = 5000 V/m. This value of conductivity is within 
the range determined in the experiments and is the 
same constant value that was used for the perfect gas 
and equilibrium air computations. 

The computed streamwise variation of averaged ve- 
locity for the nonequilibrium seeded-air calculations 
is shown in Fig. 5. The velocities are averaged across 
the channel cross section and normalized using the 
entrance velocity to be consistent with the experi- 
ment. In the experiment, the velocities were obtained 
by measuring the voltage generated by the flow at 
the last electrode pair (19) which is unpowered. This 
procedure inherently involves an averaging of the ve- 
locity profile. The numerical results indicate an in- 
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crease in the averaged velocity of about 27% with 
E y = 5000 V/m and this agrees exactly with the ex- 
perimental value of 27%. 

The computed streamwise variation of averaged ve- 
locity for the different chemistry models is shown in 
Fig. 6 for E y = 5000 V/m. Both the equilibrium air 
and the nonequilibrium seeded-air chemistry mod- 
els give similar results. The centerline variation of 
static temperature for the nonequilibrium seeded-air 
model is given in Fig. 7 for the different electric field 
strengths 


Test Case 2: NASA Ames MHD Run 16 

C Vcap. . = 380V) 

The dimensional flow parameters for this test case 
are: 

Po = 9.92 x 10 5 N/m 2 
T 0 = 5560 K 
T w = 300 K 

cr e = 130 mho/m (or variable) 

B z = 0.0, 0.92 T 

E y = 0, 4309, 5300, 6000 V/m 

This test case was also computed using several dif- 
ferent electric field strengths in order to properly 
simulate the experiment. The electric field strength 
of 4309 V/m corresponds to the voltage drop mea- 
sured in the central inviscid core flow. The com- 
puted streamwise variation of static pressure for the 
nonequilibrium seeded-air calculations is shown in 
Fig. 8 for the different electric field strengths. The ex- 
perimental pressure variation agrees with the numeri- 
cal result with an electric field strength of 5300 V/m. 
The computed streamwise variation of static pressure 
for the different chemistry models is shown in Fig. 9 
for E y = 6000 V/m. The different chemistry models 
produce similar results. 

The averaged electrical conductivity at the center 
of the powered portion of the channel (electrode pair 
10) was found to be 142 mho/m for E y = 5300 V/m. 
This value is at the upper end of the range deter- 
mined in the experiments and is higher than the con- 
stant value of 130 mho/m used in the perfect gas and 
equilibrium air computations. 

The computed streamwise variation of averaged ve- 
locity for the nonequilibrium seeded-air calculations 
is shown in Fig. 10. The numerical results indicate 
an increase in the averaged velocity of about 38% for 
E v = 6000 V/m and this agrees closely with the ex- 
perimental value of 39%. The computed streamwise 


variation of averaged velocity for the different chem- 
istry models is shown in Fig. 11 for E y = 6000 V/m. 
Once again, the equilibrium air and the nonequilib- 
rium seeded- air models give similar results. The cen- 
terline variation of static temperature for the equi- 
librium model is shown in Fig. 12 for the different 
electric field strengths. 

Concluding Remarks 

In this study, a new 3-D parabolized Navier-Stokes 
algorithm with nonequilibrium seeded-air capability 
has been developed to efficiently compute MHD flows 
in the low magnetic Reynolds number regime. The 
new algorithm has been used to compute the flow 
in the NASA- Ames experimental MHD channel for 
Runs 15 and 16. The numerical results are in good 
agreement with the experimental results. 
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Figure 1: Cross section of NASA Ames MHD Channel 



Figure 2: Schematic of powered portion of MHD channel 
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Figure 4: Streamwise variation of static pressure for Run 15 ( E y = 5000 V/m) 
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Figure 6: Streamwise variation of averaged velocity for Run 15 (E y = 5000 V/m) 


12 

American Institute of Aeronautics and Astronautics 





Figure 8: Streamwise variation of static pressure for Run 16 (seeded-air model) 
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Figure 10: Streamwise variation of averaged velocity for Run 16 (seeded-air model) 
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Figure 12: Streamwise variation of static temperature for Run 16 (equilibrium air model) 
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